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The density matrix renormalization group (DMRG) method introduced by White for the study of 
strongly interacting electron systems is reviewed; the method is variational and considers a system 
of localized electrons as the union of two adjacent fragments A,B. A density matrix p is introduced, 
whose eigenvectors corresponding to the largest eigenvalues are the most significant, the most prob- 
able states of A in the presence of B; these states are retained, while states corresponding to small 
eigenvalues of p are neglected. It is conjectured that the decreasing behaviour of the eigenvalues is 
gaussian. The DMRG method is tested on the Pariser-Parr-Pople Hamiltonian of a cyclic polyene 
(C-ff)jv up to N = 34. A Hilbert space of dimension 5. x 10 18 is explored. The ground state energy 
is 10 _3 eV within the full CI value in the case N = 18. The DMRG method compares favourably 
also with coupled cluster approximations. The unrestricted Hartree-Fock solution (which presents 
spin density waves) is briefly reviewed, and a comparison is made with the DMRG energy values. 
Finally, the spin-spin and density-density correlation functions are computed; the results suggest 
that the antiferromagnetic order of the exact solution does not extend up to large distances but 
exists locally. No charge density waves are present. 



I. INTRODUCTION 

A few years ago White introduced in the study 
of electron correlation a new and powerful numeri- 
cal method: the density matrix renormalization group 
(DMRG). The method provided extremely accurate re- 
sults in the case of the one-dimensional Heisenberg and 
Hubbard models aH, Hubbard-like models with bond 
alternations and recenjtly has been applied to some two 
dimensional models Lru. 

DMRG is a new variational method that promises to be 
very useful in quantum chemistry. It deals with the main 
difficulty of this kind of calculations, i.e. the exponen- 
tial increase of the dimension of the Hilbert space with 
the size of the system, in a new, direct and efficient way. 
While the usual packages of ab initio quantum chemistry 
cut the dimension of the Hilbert space by neglecting the 
coefficients of the configuration interaction expansion be- 
low a certain threshold, the DMRG obtains an analogous 
result with a different strategy. A system of localized 
electrons is partitioned in two blocks A,B (sometimes B 
is called "environment" or "universe") and all the many- 
electron states corresponding to situations in which the 
population of the two blocks is unphysical (e.g. all elec- 
trons in A, no electrons in B) are automatically truncated 
by the formalism. A density matrix is introduced, whose 
eigenvectors, corresponding to the larger eigenvalues, are 
the most significant, the most probable states of A in the 
presence of B. 

These states are retained, and states corresponding to 
very small eigenvalues are neglected. The two blocks are 



taken initially small and increase their size in the course 
of the calculation. As a result of the systematic trunca- 
tion mentioned above, the time of computation does not 
grow more than the fourth power of the size of the system, 
keeping constant the number m of-retained eigenvectors 
of the density matrix. The errora is an exponentially 
decreasing function of m. 

The method is especially suited to treat systems with 
translational or reflection invariance, since in an interme- 
diate stage of the calculation wave functions suitable to 
describe the block B can be obtained simply by transla- 
tion (or reflection) from those of block A. 

A good candidate in order to test the method in 
quantum chemistry is provided by the Pariser-Parr-Pople 
model of conjugated polyenes. Many considerations are 
in favour of this choice: 

• A cyclic polyene (CH)n with the carbon atoms at 
the vertices of a regular polygon is "translation- 
ally" invariant (here translation means a rotation 
of the circle circumscribed to the polygon); hence 
the simplification mentioned above can be applied. 

• Exact full. c o nf iguration interaction calculations are 
available MLil, and we can compare the DMRG 
ground state energy values with these results. The 
comparison can be made up to N — 18. A further 
comparison can-be made with the coupled cluster 
(CC) method B. However, the DMRG method 
is much more powerful; we have computed with- 
out much effort ground state energy values up to 
N = 34 carbon atoms. The full CI Hilbert space 
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corresponding to N = 34 has dimension equal to 



34V 
17) 



5.44 x 10 18 . 



• Trans-polyacetylene presents interesting experi- 
mental and theoretical problems: the bond alterna- 
tion (and in particular the values of the two bond 
lengths) can in principle be deduced by ab initio 
computations but this problem meets considerable 
difficulties. Recently an interesting approach to 
the problem of the dimerization of polyacetylene 
using the DMRG method has been put forward 
by M.B.Lepetit and G.M. Pastor L3; these last au- 
thors treat accurately the hopping term allowing a 
dependence on the distance between the C atoms 
and describe the electron interaction by a Hubbard 
term. Therefore it would be of interest to extend 
their work by substituting a PPP interaction to the 
Hubbard interaction. In the present paper we show 
that this extension is possible (but we do not derive 
the hopping term from ab initio calculations). 

• The unrestricted Hartree-Fock solutions of the PPP 
model Hamiltonian p*^seHt spin density waves and 
charge density waves EJii3. It is of interest to know 
whether or not these waves persist after a more pre- 
cise variational approximation to the ground state 
(like the DMRG) is performed. 

The paper is organized as follows: in Sec. 2 the PPP 
Hamiltonian is written down and the DMRG method is 
reviewed. In particular we point out some mathemat- 
ical aspects of the DMRG method that usually are not 
sufficiently emphasized. In Sec. 3 the properties of the un- 
restricted (spin density wave) Hartree-Fock solution are 
briefly discussed. In Sec. 4 the numerical results and the 
conclusions are presented. 



II. THE PPP HAMILTONIAN AND THE DMRG 
METHOD. 

The Pariser-Parr-Pople Hamiltonian of the it electronic 
1 of a cyclic polyene CnHn can be written as 



H = /3 
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where E^ v are the generators of the unitary group 
summed over spin, and n M = is the occupation num- 
ber of the site /i; /3, 7^ are parameters of the model, and 
< u,v > denotes summation restricted to nearest neigh- 
bor. We limit ourselves to the series N — 2n = 4^+2, v = 
1,2,..., where N denotes the total number of electrons 
which is equal to the total number of sites. According 
to ref. [ |l^| we take j3 = —2.5 eV, and for the Coulomb 
repulsion we use the Mataga-Nishimoto prescription E3: 
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where 70 = 10.84 eV, d M „ denotes the distance between 
the vertex /1 and the vertex v of a regular polygon of N 
sites and is given by 



sin(/i-^)f 



(2.3) 



and b, the nearest- neighbor separation, is equal to 1.4 A. 

Let's see now how the DMRG method can be applied to 
the PPP model. We will also-ceview briefly the principal 
formal and physical ideasEroEj that are behind this new 
and powerful numerical method. 

Let A and B denote two adjacent subsets of respec- 
tively Na, Nb sites. The method consists of two parts: 
step 1, called the "infinite system method" and step 2, 
called the "finite system method". In step 1, Na + Nb < 
N, Na = Nb and Na, Nb are progressively increased 
up to reach the condition Na + Nb — N, while in step 
2, we have always Na + Nb = N, with variable Na 
and Nb- For instance in step 1 we can have N = 18, 
A = {1,2,. .6}, B = {7, 8, ..12}, in step 2 we can have 
A = {1,2,3,4}, B = {5,6, ....18}. The main task of 
the method is to find a reduced set of "localized" many 
particle states for subsets (blocks) A and B suitable to 
describe the union A (J B. 

Let us denote by A + , B + polynomials in the creation 
operators corresponding to sites in A, B, respectively. 
Let |0) denote the vacuum, and let \a) — ^4 + |0), \b) = 
B + \0). Clearly |o), |6) represents states of electrons local- 
ized in different subsets. We can form the state A + B + 10) ; 
this state is similar but not identical to the tensor product 
I a) ® \b) since the operators A + , B + do not necessarily 
commute. We use the notation \a)\b) to denote the com- 
pound state A + B + \Q). Clearly, varying the polinomials 
A + , B + in all possible independent ways, the states \a)\b) 
generate the whole Hilbert space. 



TABLE I. Energy results: the enerj 
via Restricted and Unrestricted HF. 



jies (in eV) calculated 
FCP and DMRG are 



compared for different values of N. indicates the number 
of states kept in block A during the the n-th DMRG iteration. 



N Erhf 



Euhf 



Efci 



Edmrg m 



(1.2) 



,0) 



6 
10 
14 
18 
22 
26 
30 
34 



-11.358325 
47.441467 
-23.731302 
-30.101389 
-36.513220 
-42.950070 
-49.403281 
-55.867856 



-11.358325 
■17.910422 
-24.924267 
■32.007998 
■39.105943 
■46.207715 
■53.310920 
■60.414852 



-12.722033 
-20.060503 
-27.671391 
-35.385430 



-12.722032 
-20.060503 
-27.671333 
•35.384861 
-43.145027 
■50.928028 
-58.715323 
-66.509902 



256 
256 
256 
256 
256 
256 
200 
200 



512 
512 
512 
512 
512 
512 
400 
400 



a From Ref. [ MOM 
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Suppose that we have found an exact or approximate 
ground state \ip) of Na+Nb electrons in the subset A [J B 
of the chain; let us expand as: 

|V> = 5>z./A+B+|0) (2.4) 
u 

where {A^~|0)} denote a complete orthonormal set of 
states of electrons localized in A, and {£?j|0)} is 
an analogous complete set of states corresponding to 
B. For instance, initially we can have A^Bj |0) = 
a^a^...a^ A 6~^6^...6^ s |0) where the a+ , b + create elec- 
trons in A, B respectively ; in this case the numbers ipu 
are the usual configuration interaction (CI) coefficients. 
In principle the sums Yin Ylj run over 4 > 4 s states 
respectively, since the occupation numbers nj , of a site 
can have four possible values: (0, 0), (1, 0), (0, 1), (1, 1). 
However the number of spin up electrons and the num- 
ber of spin down electrons are good quantum numbers 
and can be fixed; we can choose states Af |0), Bj |0) with 
fixed numbers of spin up and spin down electrons, and 
the coefficients tf>u vanish unless this conservation law 
is fulfilled. Furthermore, during the iteration procedure, 
the number of states will be truncated; therefore in the 
expansion (2.4) we keep in general only tua states for the 
block A and tub states for the block B. In the following 
we shall assume that the coefficients ipu are real. 

The main mathematical tool of the DMRG theory is 
provided by the following density matrix: 

m B 

Pip = £ tirfi'J = {^ T )n> ( 2 - 5 ) 
j 

The dimension of the matrix p is vtla x to^; however, 
because of the number conservation laws described above, 
the matrix is actually in block form: the number of up 
and down electrons of the states I and I' must be the 
same. Furthermore p is a non negative square matrix. 

Let us first make some simplifying assumptions, that 
will be relaxed in the following. Let's assume that the 
blocks A and B are described by the same number of 
states (rriA = hi-b), so that the matrix ip is a square 
matrix. Denoting by S the square root of p (p = i/ji[> t = 
S 2 , S = ps), we have the polar decomposition 

ip = SUi (2.6) 

where U± is an orthogonal matrix. We diagonalize S by 
writing S = UDU T , where U is an orthogonal matrix 
and D is diagonal. Therefore we can write 

^ = UDU T Ui = UDV T (2.7) 

where V is an orthogonal matrix, and p = UD 2 U T . 

Actually formula (2.7) holds for any rectangular matrix 
niA x nriB ip (see e.g [ p2[). U and D are square matrices 
niA x tua and V T is tua x tub- These matrices verify 
the conditions: 



U T U = I ; VV T = I ; 

^ T = UD 2 U T ; V 7 V = VD 2 V T (2.8) 

and D is diagonal and non-negative. Let us denote by 
D a the eigenvalues of D. Substituting (2.7) into (2.4) we 
obtain: 

m A 

|V>> =^D a \u a )\v a ) (2.9) 

a 

where 

rriA TUB 

K) = J2 U Ia A+\0) , \v a ) = VjaB$\0) (2.10) 
i j 

What is the meaning of \u a ), \v a )7 They represent states 
of the subsystems A, B, such that the probability for the 
whole system A\J B to be found in the state |ti a )|u a ) is 
D 2 a . The main idea of the DMRG method consists in ne- 
glecting, in Eq.(2.9), all eigenvalues D a below a certain 
threshold which amounts to keeping only a small number 
til of terms in the sum (2.9) and using the corresponding 
states \u a ) as a basis for the description of block A. Since 
Tripij} T — Y^2 A = I; this approximation is good if the 
probabilities D 2 a have a sufficiently rapid decrease to zero, 
so that Xm=i — !• At the best of our knowledge, all 
numerical experiments performed so far (see, e.g. ref. 
[ ^3|j2l|l ) confirm this rapid decrease of the probabilities 
D a . Let's give an heuristic argument for this decreasing 
behaviour. Suppose that A^ creates N A electrons in the 
Na sites of the block A, and Bj creates N B electrons 
in the Nb sites of the block B. In absence of the inter- 
action, usual statistical mechanics arguments prove that 
the probabilities D 2 a are strongly peaked about the pop- 
ulations N A = Na, N b = Nb (which correspond to a 
density of one electron per site); this is analogous to the 
classical result in statistical mechanics stating that the 
probability of distributing a large number of molecules 
in two communicating volumes is strongly peaked about 
a distribution with equal density in the two volumes. 

TABLE II. Correlation energies: The correlation energies 
per electron (in eV) of the FCI a and DMRG solutions with 
respect to the Restricted and Unrestricted HF approximations 
are compared for different values of N. 

E-E(RHF) E-E(UHF) 



N N 



N 


FCP 


DMRG 


FCP 


DMRG 


(1.2) 
m A 


(3) 


6 


-0.227285 


-0.227285 


-0.227285 


-0.227285 


256 


512 


10 


-0.261904 


-0.261904 


-0.215008 


-0.215008 


256 


512 


14 


-0.281435 


-0.281431 


-0.196223 


-0.196219 


256 


512 


18 


-0.293558 


-0.293526 


-0.187635 


-0.187603 


256 


512 


22 




-0.301446 




-0.183595 


256 


512 


26 




-0.306844 




-0.181551 


256 


512 


30 




-0.310401 




-0.180147 


200 


400 


34 




-0.313001 




-0.179266 


200 


400 



a FromRef. [|l^0. 
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Because of the central limit theorem, the peak is gaus- 
sian in the classical case; we make the conjecture that 
even in the quantum interacting case that we are con- 
sidering, this gaussian behaviour still holds, at least for 
translationally invariant systems, like the PPP model. Of 
course if the conjecture is true, it explains the strongly 
decreasing behaviour of the probabilities D\ mentioned 
above. 

Let's now proceed with the description of the DMRG 
method. Once we have a pretty good basis of ttla states 
\u a ) that describe the block A, and m B states \v a ) that 
describe the block B, the next task consists in the en- 
largement of the blocks. In the first part of the algorithm 
(infinite system method) , since Na — Nb and the system 
is translationally invariant, the states \v a ) can be simply 
obtained by translating the states \u a ). Hence we can 
concentrate our attention on the block A. 

The simplest way of enlarging the block A consists in 
adding .a, site s to A, obtaining a new block A' — A [J s. 
White u denotes by A • this new block. Denoting by 
|*i> = |0>, |» 2 ) = <4 T |0>, | S3 ) = oJjO), |« 4 > = a^oJjO), 
the states describing the site s, we have Am a vectors 

A' I + \0) = \u ai )\s 7i ), /=(a I ,7 J ) 

a 1 =l,...m A , 7 7 = 1,...4 (2.11) 

in order to describe A' = A • = A (J s. 

At the same time, we add an analogous site t to 
the block B, and we consider the vectors \vp)\ts) {(3 — 
1,2, ...m^, 5 = 1,2,3,4) in order to describe the block 
B' = B • = B[Jt. With such a basis we can now pro- 
ceed to compute the expansion (2.4) for the wavefunction 
for the new superblock A'[jB'. 

Let us use the term "local" to denote operators 
at , , referring to one site \i only, and the term " in- 
ternal to block A" to denote operators whose site indices 
belong to the block A. 

The idea is now to compute a new "effective" Hamilto- 
nian matrix H', by using the truncated basis consisting 
of the 16myimB vectors |it a )|s 7 )|^/3)|^). 

TABLE III. Comparison of approximate solutions: The 
correlation energy E ~ E ^ tHF (hi eV) of the DMRG solution is 
compared with the partial cluster analysis (\e D ) \RHF)) a , the 
Approximate Coupled Pair theory with Quadruples (ACPQ) b 
and the Approximate Coupled Pair theory with Triples and 
Quadruples (ACPTQ) b . 



N 


\e D )\RHF) a ACPQ b 


ACPTQ b 


DMRG 


(1.2) 


(3) 


6 


-0.224196 


-0.2238 


-0.2253 


-0.227285 


256 


512 


10 


-0.248723 


-0.2515 


-0.2577 


-0.261904 


256 


512 


14 


-0.256777 


-0.2649 


-0.2762 


-0.281431 


256 


512 


18 




-0.2720 


-0.2887 


-0.293526 


256 


512 


22 




-0.2763 


-0.2994 


-0.301446 


256 


512 



a FromRef. [H. 
b From Ref. [jig. 



Clearly it is easy to compute terms of the Hamiltonian 
containing local operators referring only to one of the four 
blocks A, s, B, t; these terms are known from previous 
steps of the iteration. A little more care is needed in 
order to compute terms like a^a v or n^n^ with fi, v 
belonging to different blocks (e.g. fi € A, v € s, etc.). 
For this purpose it is necessary to keep in the computer 
memory all the matrix elements of the local operators 
(u ai \al\u a2 ), (vpArivlvfa). 

The entire procedure can now be repeated: we look for 
the ground state vector ip 1 of the truncated Hamiltonian 
H', by using Lanczos's or Davidson's algorithm. A new 
density matrix ip'ip' T and new state vectors \u' a ), that 
represent states of A', are computed according to the 
analogous of the first of formulas (2.10) which now reads: 

4m_A 

K> = £tM + |o> 

4mA 

= E C/ ^K>K> <* = l,-m A > (2-12) 
i=i 

Again we do not keep all the vectors: tua 1 is generally 
less than Am a and often one puts vcia 1 = mA-, although 
this choice is not necessary. The corresponding \v'g) that 
describe B' are obtained from the \u' a ) by translation. 

In this new truncated basis we compute the matrix 
elements of all the local and internal operators relative 
to block A' and we keep them in the computer memory, 
in order to use them in next steps of the method. If, 
for example, we have an operator O internal to block A, 
it is also internal to the new block A' and we have the 
following rule to update its matrix elements: 

4mA 

J,/'=l 

4mA 

= E u 'i ai { u ai \0\u ai ,){s li \s li ,)U' I , a2 
i,i'=i 

4mA 

= E U 'l a A U °' I \ \ U ^ ll ) 5 l I l I , U 'l' a ^ 
I,I' = 1 

for ai,ot2 = 1) ■■■Tn A > (2-13) 

Two more sites are added to the blocks A' , B' , giving 
rise to new blocks A" = A' ; B" = B' •, etc. By the sys- 
tematic procedure of adding two more sites, truncating 
the basis and updating the Hamiltonian matrix at each 
iteration, systems of large size can be handled. 

A comment is in order about the choice of the two sites 
that are added and their position with respect to blocks 
A and B. We can form the superblock A • B • or the 
superblock A • • B. White suggests that the enlarged 
configuration A • B • is to be preferred to A • • B in the 
case of periodic boundary conditions, the opposite holds 
in the case of open boundary condition. A *B •, 
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N = 34 

0.15 ■ 
0.1- 



T 0.05 
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FIG. 1. Spin-spin correlation function versus the distance 
between sites. 



In fact the blocks A and B are separated by the site t 
in the case while they become adjacent by periodicity in 
A • • B. The kinetic part of the Hamiltonian (2.1) "con- 
nects" two blocks only by its border sites with operators 
aL a„, whose matrix elements are known. These ma- 
trices are "big" for blocks A and B, and "little" for the 
1-site blocks s and t, so the matrix elements of the hamil- 
tonian H' are simpler when a "big" block is surrounded 
by 1-site blocks. 

The "infinite system algorithm" is stopped when the 
number of sites of A [j B reaches the total number N of 
sites. In order to improve the accuracy of the method, 
White himself proposed a second algorithm, that we will 
briefly describe. This second algorithm takes place after 
the infinite system algorithm reaches the end. 

In the finite system algorithm, to an increase of A 
by one site it corresponds a decrease of the "universe" 
B by one site. Denoting by A x , B y blocks A and B 
with x, y sites respectively, we start with the system 
An__ 1 • Bn__ 1 • and we want to construct the systems 
An • Bn o •■> An , i • Bn , • , etc. Therefore in order 
to use the translational invariance, we need to keep in 
the computer memory all the relevant matrix elements 
of An. _ 2 , An_ _ 3 , etc., in order to be able to use the sym- 
metry and produce the matrix elements of Bn__ 2 , Bn__ 3 , 
etc. It should be noticed that when ms' < ttla 1 (this 
certainly happens when B becomes small) the rows of i/j' 
cannot be linearly independent. As a consequence, ip'ip' T 
has many eigenvalues equal to zero. 

From Eq.(2.8) we see that the vtia 1 x tua' matrix ip'tp' T 
and the smaller ros» x tub 1 matrix ip' T ?p' have the same 
non vanishing eigenvalues. In practice it is sufficient to 
diagonalize only the smallest of the two density matri- 
ces. The procedure stops when we reach the system 



Aat_3» Bi», i.e. when the block B has reduced to a 
single site. We can now increase B and decrease A; the 
subsystems A, B behave like if they were separated by a 
moving zipper. At every step we increase the accuracy of 
the states \u a ) that describe the blocks A x and after few 
oscillations of the zipper all the blocks A' x , 2 < x < N—2, 
accurately represent parts of a complete system of N 
sites, the remaining "universe" being the corresponding 
B' N _ X block. 

During this procedure, all the relevant matrix elements 
of the local operators must be stored and updated. A 
more detailed explanation ofJjhis point can be found in 
the original paper by White □. Usually one stops when 
A and B have the same length. 

III. THE UNRESTRICTED HARTREE-FOCK 
(UHF) SOLUTION. 

Let us still denote by the creation operator of an 
electron of spin a on the site /i. The creation operator of 
an electron in a symmetric Bloch orbital is given by (we 
use letters k, fci, ki.. to denote the symmetric orbitals): 

N-l 

°t = TP E e<Wfc/i < fc = 0,l,..JV-l (3.1) 

where u> = In terms of these operators, the Hamil- 
tonian can be written as 

H = 2(3 ^2 cos(wfe)a2 CT afc CT - E 

+ \ E K( - k ) a tla a t 3 + k,T a k 3 .ra kl +k,a (3.2) 
k\k 3 kar 

where all k indices run from to N— 1, the constant term 
Eq = X) 7 <i/7A" y nas been added to the Hamiltonian and 
represents the internuclear repulsion energy, and K{k) is 
given by: 

N-l 

^) = ^Ev'^ fc = 0,l,..JV-l (3.3) 

A* 

Due to the discrete rotational symmetry of the polygon, 
all indices can be taken modulo N . It is convenient to 
represent the k indices on a circle (see fig. 1 of ref . [ [l6| ) . 
The restricted Hartree-Fock orbitals are determined by 
the condition: 

N -v <k<N + v (3.4) 

which characterizes the Fermi sea F. The restricted 
Hactree-Fock (RHF)single particle energies are given 
byta 

e k = 2/3cos(wfc) + NK(0) - ^ K(k - ki) (3.5) 

feieF 
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and the total RHF energy is: 



where 



r(e(fc) - e(k + n)), and the UHF orbital 



Erhf = -E + ^[2/?cos(wfc) + e k ] 



(3.6) 



k£F 



It is known since long time that it is possible to lower 
the RHF ground state energy by considering molecular 
orbitals that are linear combinations of the orbitals <p k 
and 4>k+n corresponding to two cndpoints of a diameter 
of the circle of Fig.l in ref. [ 16 1. 

Furthermore, taking into account also the spin indices 
of the two orbitals, there are many different possibilities 
that give rise to local minima of the UHF energy. All 
these possibilities-have been carefully studied many years 
ago by FukutometZl, Paldus and CizekEj, and others, and 
give rise to charge density waves (CDW) and spin den- 
sity waves (SDW). However in the case of the Mataga- 
Nishimoto prescription for the two center Coulomb repul- 
sion integral and with the values of the parameter given 
in Sec.l, we have checked that the lowest UHF energy 
is obtained by the following BCS-Bogoliubov canonical 
transformation (which corresponds to the (A t +B t ) + case 



of ref. [|| 



7fcT = u fe fl feT + Vkdk+nl 



1H = ~u k a k[ + v k a k+nl (3.7) 

where u\ + «| = 1 , u k+n = u k , v k+n = -v k . The 
operators 7^ create UHF orbitals, since the linear com- 
bination depend on a. 

The first-order density matrix (in the pseudomomenta 
representation) is given by: 

(3.8) 

where f^(k) — u 2 , f^ 2 \k) = u k v k for k e F, and 
/( 1 )(/c) = v%, / (2) (fc) = -u k v k for k g F. In the original 
atomic-orbital basis we have the interesting formula: 



(3.9) 



where 5 = -h Yl k =o l u fe w fe|- This formula shows the exis- 
tence of SDW of antiferromagnetic type; the occupation 
numbers and (n^) are different from each other on 
the same site (when one of the two is larger than i the 
other is smaller than i) giving rise to a decrease of the on- 
site Coulomb repulsion. Furthermore no CDW appear, 
since (n^) + (n^) = 1. 

The expectation value of the Hamiltonian (2.5) can be 
easily computed by using Wick's theorem; minimization 
of (H) with respect to the coefficients u/^irgives rise to 
the following well known set of equationalaO of the BCS 
type: 



= ; 1 



la 




energies are given by 
e(k) = 2/3cos{ujk) + NK(0) - ^ K(k - q)u 2 q 

q£F 

-^K{k- q)v 2 q , for k = 0, 1, ...N - 1 

q?F 



(3.11) 



and A(fc) must fulfil the famous "gap equation": 
A(fc)= 1 1 E q eF[K(k-q) + K(k~q + n)] 

X 7 W , forfc = 0,l,..JV-l (3.12) 

If the only solution of the gap equation is the trivial 
solution A(fc) = 0, we obtain simply the RHF ground 
state. If a non trivial solution exists, the non linear 
system of equations (3. 10), (3. 11), (3. 12) can be easily 
solved numerically by an iterative method. Starting with 
A(q) — constant and e(k) = e(fc), we solve the gap equa- 
tion (3.12) by iteration. Usually 30 — 40 iterations will 
suffice. The solution A(fc) is substituted into Eqs. (3.10), 
(3.11); in this way we obtain a set of approximate orbital 
energies ei(fc). The entire procedure is repeated substi- 
tuting in the right hand side of Eq. (3.12) the solution 
A(/c) and = i(ei(fc) — ei(k+ n)), etc., until the entire 
set of equations is fulfilled with sufficient accuracy. 
The UHF ground state energy of the model is given by 



JV-l 



E UH F = - J2 K ( k 
fc 1: fc 2 =0 

/ (1) (fcl)/ (1) (fc 2 ) + / (2) (fcl)/ (2) (fc 2 ) 

-E + - K(0) N 2 +4f3^2 cos(wfc) / (1) (fc) 



(3.13) 



The antiferromagnetic long-range order of the UHF so- 
lution appears also in the height of the peak of the mag- 
netic structure factor: 

^)=Ef=o 1 e- i ^' fe <^(j)^(0)) 



+^ f 5 2 -^Er=o i / (i) (?)/ (i) (9-fc) 

-^El-of (2) (<l-k-n)fV(q) 



which is reached for k — y. We have 



S(- 



,N s 



(3.14) 



(3.15) 
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and this quantity scales like N for large N. 



(3.10) 



G 



0.3 



0.2 



N = 34 



i 0.1 




-0 . 1 



FIG. 2. Density-density correlation function versus the 
distance between sites. 



IV. NUMERICAL RESULTS AND 
CONCLUSIONS 

In Table | we show the energy results calculated with 
the DMRG method up to N — 34, and we compare them 
with RHF, UHF and FCI energies (the FCI energies are 
available only up to N — 18). 

We see that the relative error of the DMRG solution 
with respect to FCI is only 2.1 x 1CT 6 for N = 14 and 
1.6 x 10 -5 for N = 18, which is a quite satisfactory result. 

Table [n] shows the correlation energy per electron of 
the FCI and DMRG solutions with respect to the RHF 
and UHF approximations. 

The DMRG method compares favourably with the 
Coupled Cluster method; in Table III the correlation en- 
ergies (E — Epjjp)/N are compared with coupled cluster 
results of ref. [ y>H. The DMRG energy is slightly 
lower than the Approximate Coupled Pair with Triples 
and Quadruples (ACPTQ) value. All calculations were 
performed iterating the DMRG algorithm three times 
(the first iteration uses the infinite system method, the 
second and third iterations use the finite system method). 
We stop when A and B have the same length. 

In the first iteration the size of the system grows, but 
the potential between two sites is kept equal to its fi- 
nal value, i.e. to the value attained when the number of 
sites of the polygon is N. Generally we keep 256 states in 
block A during both the first and the second iteration; in 
order to achieve a better convergence, during the third 
iteration we keep 512 states. In the heaviest calcula- 
tions (N = 30, 34) we keep only 200 - 400 states in block 
A, due to memory /disk-space limitations. It should be 
noted that the disk-memory requirement grows with the 
number of sites even if the number of retained states is 



held constant. This is due to the long range nature of 
the interaction that forces us to keep on disk a linearly 
growing number of matrices that represent the local op- 
erators. We have checked that the disk-space grows as 
Nm 2 . 

In Fig. the spin-spin correlation function S(i — j) = 
(S z (i)S z (j)) is plotted : a short range antiferromag- 
netic order is clearly present. We have computed the 
Fourier transform S(k) which of course reaches its max- 
imum value for k = f , like the UHF-SDW solution 
(see (3.14)). However, the growth is linear with N 
for the UHF-SDW solution, but scales approximately as 
0.1398 + 0.1457LogA for the DMRG solution. Therefore 
we cannot speak of long range SDW. Also the CDW are 
ruled out by the present calculation. This can be seen 
from the graph of the density-density correlation function 



R{hj) = (n{i)n(j)) - (n{i))(n(j)} 



(4.1) 



(see Fig. |). 

Concluding, the DMRG method provides a very pow- 
erful tool for the calculation of energies and properties 
of simple many electrons Hamiltonians. It gives results 
very close to full CI results and is able to handle Hilbcrt 
spaces of very large dimension. It would be of great in- 
terest to apply the method to a realistic many electrons 
Hamiltonian, possibly after a previous localization of the 
occupied and virtual orbitals. However, this program 
meets with some difficulty because of the large number 
of matrices that must be kept when the four orbitals of 
the interaction term belong to different blocks. 
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